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The presence of attractive interaction between fermions can lead to pairing and superfluidity in 
an optical lattice. The temperature needed to observe superfluidity is about a tenth of the tunneling 
energy in the optical lattice, and currently beyond experimental reach. However, at strong coupling 
the precursors to global superfluidity should be visible at achievable temperatures, in terms of 
fluctuating domains with strong pairing correlations. We explore this regime of the attractive two 
dimensional fermion Hubbard model, in the presence of a confining potential, using a new Monte 
Carlo technique. We capture the low temperature inhomogeneous superfluid state with its unusual 
spectral signatures but mainly focus on the experimentally accessible intermediate temperature 
state. In this regime, and for the trap center density we consider, there is a large pairing amplitude 
at the center, spatially correlated into domains extending over several lattice spacings. We map out 
the thermal evolution of the local density, the double occupancy, the pairing correlations, and the 
momentum distribution function across this phase fluctuation window. 



The development of optical lattice methods for ultra- 
cold atoms has opened a new vista in the study of cor- 
related systems [US], allowing clean controllable reali- 
sations of strongly interacting quantum lattice models. 
Experimental achievements include the realisation of the 
superfluid (SF) to Mott insulator transition [5 in the 
Bose Hubbard model, and, for repulsive fermions, the 
observation of Fermi surface [6 , and the Mott insulat- 
ing phase [7J [8]. For attractive interactions, there has 
been the evidence of superfluidity [9 , a possible FFLO 
state [10] , and anomalous expansion of the Fermi gas [IT] . 
The realisation of an antiferromagnetic state p~2j [13] in 
the repulsive Hubbard model and of superfluidity in the 
attractive model [14] is still awaited. The problem is with 
the achievable temperature [T5] . 

Techniques available to date can reduce the entropy, 
£, per particle of a Fermi gas to ~ log e 2 « 0.7. The 
associated temperature is C(t), where t is the tunnel- 
ing energy between neighbouring wells in the periodic 
potential. The observation of superfluidity in the attrac- 
tive Hubbard model will require cooling to ksT/t ~ 0.1. 
The corresponding entropy is S ~ 0.1 [15] . almost an 
order of magnitude below what is currently achievable. 
What signatures would one expect of attractive interac- 
tions at accessible temperatures? At weak coupling the 
state above T c is a normal Fermi liquid but at strong 
coupling a large pairing amplitude survives to T > T c 
and entropy levels S ~ log e 2. This is an inhomogeneous, 
thermally fluctuating, short range correlated state, with 
striking measurable properties. 

In this paper we solve the attractive ('negative L 7 "') 
Hubbard model on large two dimensional lattices in the 
presence of a confining potential. Our main results are 
the following: (i) We access the superfluid ground state, 
the thermal transition, and a wide temperature window 
over which the Fermi system has large pairing amplitude 
but no global phase coherence, (ii) We illustrate how fluc- 



tuating filamentary 'superfluid' regions survive far above 
T c , to temperatures ksT/t ~ 0(1), and leave signatures 
on the spatial patterns and spectral features. We demon- 
strate these using a new method that handles the strong 
coupling non-perturbatively and treats the spatial inho- 
mogeneity and strong thermal fluctuations exactly. 

We study the attractive Hubbard model in the pres- 
ence of a harmonic potential Vi in two dimensions: 

H = -t ^ C \a C jcr + y^(Vi ~ lA^cr ~ U ^2 n *t n 4 (!) 
(ij)cr icr i 

The first term denotes the nearest neighbour tunneling 
amplitude of fermionic atoms on the optical lattice, the 
confining potential has form Vi = Vb(xf + y?) : fi is the 
chemical potential, and U > is the strength of attrac- 
tive on-site interaction, Xi and yi are measured in units 
of lattice spacing ao- 

The spatial variation in mean value, and the ther- 
mal fluctuation about the mean, of the amplitude and 
phase of the order parameter are crucial in describing the 
physics of this system. Unbiased calculations in the ho- 
mogeneous limit employ determinant al quantum Monte 
Carlo [T5HT7] (DQMC) to access finite temperature prop- 
erties, but are typically limited to 10 x 10 lattices. That 
is inadequate to clarify the interplay of correlation effects 
and inhomogeneity. 

We use a strategy used earlier on moderately sized 
systems [HI [19], augmented now by a 'traveling clus- 
ter' (TCA) [20] Monte Carlo technique that readily al- 
lows access to system size ~ 32 x 32. We first decouple 
the Hubbard term in the pairing channel by using the 
Hubbard- Stratonovich (HS) transformation. We use the 
static HS (sHS) approximation [18], z.e, retain spatial 
fluctuations of HS fields but ignore the time dependence. 
This leads to the following effective Hamiltonian: 

Heff = h + J2( A 4A + A * c ^t) + E (2) 
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where H = -t E<^> 4* c jv + T,ia( v i lA n iv and A i = 
|A^|e 2 ^ is a complex scalar classical field. This model 
allows fluctuations in both the amplitude and phase of 
the pairing field A^, and the fermions propagate typically 
in an inhomogeneous background defined by A^. 

To obtain the ground state, and in general configura- 
tions {|A^|,^} that follow the distribution P{|A^|,^} ex 
Tr c ^e~P Ref f , we use the Metropolis algorithm to up- 
date the |A| and variables. This involves solution of 
the Bogoliubov-de Gennes (BdG) equation [21] for each 
attempted update. For equilibriation we use the traveling 
cluster algorithm [20 , diagonalising the BdG equation on 
a 8 x 8 cluster around the update site. Global proper- 
ties like pairing field correlation, quasiparticle density of 
states, e£c, are computed via solution of the BdG equa- 
tion on the full system for equilibrium configurations. 

We explored the system at U/t = 2, 6, 12 and the 
maximum (system corner) potential V c ~ 1^*2* (L/2) 2 = 
E//2, J7, 2/7, where the system size is L x L. This enables 
us to systematically study the evolution from weak to 
strong coupling, as well as weak to strong confinement. 
We focus on the strong coupling, strong inhomogeneity 
case, U = 12, V c = 24 in this paper and will discuss the 
larger parameter set separately [22] . 

In the absence of the confining potential the model is 
known p~6j H7] to have a superfluid ground state for all 
densities n / 1, while at n = 1 there is the coexistence 
of superfluid and density wave (DW) correlations. The 
SF ground state away from n = 1 evolves [23] from a 
Bardeen- Cooper- Schrieffer (BCS) state at U/t <C 1 to a 
Bose-Einstein condensed (BEC) state of 'molecular pairs' 
at U/t ^> 1. The ground state can be reasonably accessed 
within mean field theory (MFT) but the finite temper- 
ature predictions of MFT becomes increasingly inaccu- 
rate with increase in U [24 . This is due to the sep- 
aration of scales between 'pair formation' temperature, 
Tf ~ 0(U), and pair condensation, ie, superfluidity, 
which is T c ~ t 2 /U. MFT captures Tf but wrongly iden- 
tifies it with the superfluid transition. 

For U/t < 1, the state at T > T c is an uninteresting 
weakly correlated Fermi liquid. As U/t increases, T c (in 
two dimensions) peaks at U/t « 5, while Tf continues 
to grow. A wide 'non-Fermi liquid' window opens up 
between T c and Tj, and the system behaves like a (hard- 
core) Bose liquid [24] for low temperature and U/t ^> 1. 
In this U ^> t regime, increasing T leads to gradual disso- 
ciation of the 'bosons', and paired and unpaired fermions 
exist in equilibrium. 

The confining potential promotes an inhomogeneous 
density profile [25], with a peak at the trap center. If the 
average density near trap center is n = 1, it can lead to 
a local DW pattern [26] [27] while the SF would show up 
away from the center where n < 1. If the total particle 
number is sufficiently small so that even the central den- 
sity is < 1, the entire system is an inhomogeneous SF 
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FIG. 1: Colour online: (a). Comparison of the T C (U) at 
n — 0.7 on 'flat' 10 x 10 lattices, obtained by two differ- 
ent techniques, the DQMC and our static HS approximation, 
(b). Schematic phase diagram in the flat case, indicating the 
window between the SF transition (T c ) and 'pair formation' 
scale (Tf) that grows with increasing U. The y scale is loga- 
rithmic to include the widely different scales of T c and Tf . The 
vertical dotted line indicates the T dependence that we ex- 
plore, (c). The growth of superfluid correlations at U/t = 12, 
V c /t = 24, as indicated by the zero momentum component of 
the pairing field correlation. Both the onset temperature and 
T = value are suppressed at lower N. System size 32 x 32. 



at low temperature. In our strong confinement problem 
we will focus on the case where the total particle number 
N & 150. This leads to a trap center density n = 0.9, 
optimising the T c for our choice of U and V c . For weaker 
confinement a larger N can be used. 

Fig.l.(a) compares the result of full DQMC calculation 
[T5] with that of the sHS scheme implemented via TCA. 
The comparison on 10 x 10 lattices, with Vo = and 
n ~ 0.7, shows that our method captures the overall 
trend in T C (U) and is even quantitatively accurate. In 
our understanding the agreement is due to the inclusion 
of the key thermal phase fluctuations within the static 
HS theory. This gives us confidence in the method when 
applied to large lattices and the presence of a potential. 
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Fig.l.(b) highlights the two key scales from the Vb = 
problem, at n = 0.7, that are relevant for us: (i) the non- 
monotonic T c scale whose maximum is roughly 0.1 St and, 
(ii) the pair formation temperature Tf as defined below. 
We know that the pair binding energy is 0(U) when U/t 
is large. To fix the prefactor, we define Tf = 2A p (0)/3.5, 
where 2A^(0) is the T = gap in the quasiparticle spec- 
trum. For U/t <C 1, both superfluidity and the pairing 
amplitude vanish when T = 2A^(0)/3.5 and T c = Tf by 
definition. At strong coupling 2A^ ~ J7, so Tf oc U. 

Fig.l.(c) shows the growth in the Q = {0,0} 
component of the pairing field correlations, D(Q) = 
|A i ||A j |cos(0 i -^)e iQ ' (ri " rj ' ) . This is non zero when 
the amplitude |A^| is finite over some region and the 
phases 6^ are correlated. This in turn promotes a non 
zero value of ((cj^cj^)) and a finite value for %(ij,T) = 

(( c lt c l| c ii c it))- We therefore use D({0,0},T) as indica- 
tor of the SF transition. This is shown for choices of par- 
ticle number N that lead to trap center density rii < 1. 
The T c scale, at N ~ 150 is roughly 0.12t. At lower 
N the onset temperatures are lower and the strength of 
pairing field correlation at T = is also smaller. 

Accessible temperatures are still ^> T c max ~ 0.1 St so 
one would have to look for non trivial interaction effects 
at T > T c . At weak coupling the T > T c state is unin- 
teresting. However, for U/t > 5, the window between T c 
and Tf is wide and well accessible with present cooling 
techniques. We highlight the particularly wide window 
at U = 12 that reaches from T/t - 0.1 - 3. 

We chose \i such that the maximum density, which 
occurs at the trap center, was always less than 1, and 
the ground state of the system does not involve any DW 
order. Our ground state is always an inhomogeneous 
SF. Fig. 2, left column, shows the spatial patterns in the 
'ground state' (T = O.OOlt). For our choice of /i, the 
density at the center is ~ 0.9. The double occupancy 
di = ((n^n^)) follows a similar profile and is almost 
double the noninteracting value (n^/2) 2 due to the strong 
interaction. The pairing field amplitude is also largest at 
the center (in the uniform case the pairing amplitude in- 
creases with n from n = to n < 1). The phase correla- 
tions are near perfect in the ground state, as the bottom 
row, left column indicates. 

The central and right column in Fig. 2 highlights the 
thermal evolution, with the results averaged over 40 con- 
figurations. The middle column is for T/t = 0.09, and 
the right for T/t = 0.50. The global order parameter 
for superfluidity vanishes at T/t ~ 0.12 but, as expected 
at large U/t, the pairing amplitude still survives. From 
T = 0.001* (left column) to T - T c (middle column) 
there is no significant change in the density pattern, the 
double occupancy, or the pairing amplitude. The pairing 
correlation (bottom row) is still dominantly positive at 
T/t — 0.09 but with hints of small (minority) domains. 

At the highest temperature, T = 0.50t, right column, 




FIG. 2: Colour online: Temperature dependence of spatial 
patterns. The results are all for N « 150. The first column 
is for T/t — 0.001 (the ground state), the second column for 
T/t = 0.090 (roughly below T c ) and the third column at T/t = 
0.50, deep in the fluctuating regime. The first row shows rii, 
the second shows the double occupancy di, the third shows 
the pairing field magnitude |A;|, the fourth shows the nearest 
neighbour pairing field correlation. System size 32 x 32. 

where we expect the average entropy to be ~ 0.5 per 
particle, the rii pattern is significantly broader and the 
associated di is more diffuse (with a slightly lower trap 
center value). The pairing amplitude is still significant, 
although the averaging has not completely restored the 
circular symmetry. The pairing correlation reveals a 
strong short range feature, and a filamentary pattern 
with lengthscale ~ 5ao- The correlations are stronger, 
overall, near the central part, but now have some strength 
towards the periphery also due to the density broadening. 

A direct measure of the correlated phase is the quasi- 
particle density of states (DOS), N(u) = ((^Z n S(uj — 
E n ))), shown in Fig.3 for T/t = 0.001, 0.25, 0.50. E n are 
the BdG eigenvalues in the equilibrium A^ background. 
At T = 0.001 there are three noteworthy features: (i) the 
pairing gap ~ U , (ii) the coherence peaks at the gap 
edge, reminiscent of 'flat' systems, and (hi) the 'spiky' 
features that arise from the quantisation of energy levels 
in this 'stiff' trap. We have checked that the sharp lev- 
els survive, but become more regularly spaced, even in 
the non-interacting case. It is now possible to measure 
the spectral function via photoemission [28] in cold atom 
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FIG. 3: Colour online: Quasiparticle density of states for 
T/t = 0.001, 0.25, 0.50. T c « 0.12, so the results are approx- 
imately for T = 0, 2.5T C and 5T C . The gap and the spiky fea- 
tures survive over this temperature window and are related, 
respectively, to 'preformed pairs' and the level quantisation 
due to the confining potential. System size 32 x 32. 

experiments. 

We observed that with increasing T the coherence fea- 
tures get wiped out and vanish by the time T ~ T c /2. 
The pairing gap, however, survives but with two modi- 
fications. The region over which N(uj) = now shrinks 
(the gap lessens) but with the loss of the coherence peaks 
there is a loss in band edge spectral weight. The quan- 
tised features still survive but are more diffuse. This 
gapped spectrum is visible even at T/t = 0.50, like the 
Mott gap in the positive U Hubbard model. 

Let us discuss the momentum distribution function 
n{k x ,k y ), Fig. 4, as the final signature of correlation 
physics. This can be measured from the velocity dis- 
tribution of the gas by switching off the trap. The left 
panel in Fig. 4 shows n(k x , k y ) for the ground state of the 




FIG. 4: Colour online: Momentum distribution function, 
n(k x ,k y ). The left panel shows the n(k x ,k y ) for a non- 
interacting fermion gas in the trap, with N = 150 and T = 0. 
Although there is expectedly no sharp 'Fermi surface', due 
to the background potential, n(k x , k y ) has strong momentum 
dependence. Middle panel is for trapped interacting fermions 
at T/t = 0.001. Right panel is the trapped interacting system 
at T/t = 0.50 



non- interacting trapped gas (at same TV). While there 
is no Fermi surface (FS) there is a strong momentum 
dependence. n{k x ,k y ) is very distinct in the interacting 
case: flat and broad with only a weak central peak in the 
ground state, and essentially flat at T/t = 0.50. Tuning 
the Feshbach resonance across the BCS-BEC crossover 
should observe this broadening. 

Conclusions: We have studied the attractive Hubbard 
model at strong coupling in the presence of a harmonic 
confining potential. Our non perturbative results high- 
light the destruction of global superfluid order at fairly 
low temperature but the survival of nanoscale fluctuating 
'paired' regions to high temperature. They leave an im- 
print on the spectral density and momentum distribution 
which serve as precursors to global coherence. 
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